Using ggplot faceting

Please follow the instructions in the README file before starting this tutorial.

---

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

---

Loading the packages

library(tidyverse)

---

Plotting with ggplot

Plotting genome coverage

• Erythromycin resistance can be conferred by 2473 bp plasmid encoding repL and ermC
• First lets load the overall genome coverage statistics for each isolate

coverages.NC_017763 <- read_tsv('data/stats.coverage.NC_017763.txt', show_col_types = FALSE)

• First lets load the read depths for each isolate mapped to the reference plasmid

depths <- 'data/NC_018969.tsv'

data <- read_tsv(file=depths, show_col_types = FALSE)
data

• This data file contains the number of reads mapping at each genome position (POS) of the reference (NC_018969.1). The columns are labelled by the name of the alignment file (*.bam)
• Lets rename the first column, as it begins with a # which can cause problems when specifying columns

data %>%
  rename(name='#CHROM')

• Lets also pivit the data to produce the long form of the table

data %>%
  rename(name='#CHROM') %>%
  pivot_longer(cols=c(-name, -POS), names_to = 'isolate', values_to = 'depth')

• And we need to extract out the isolate name from the alignment file paths, we can do this with a combination of te mutate() and str_replace() functions
• we can also save the data into a variable called long_data

long_data <- data %>%
  rename(name='#CHROM') %>%
  pivot_longer(cols=c(-name, -POS), names_to = 'isolate', values_to = 'depth') %>%
  mutate(isolate=str_replace(isolate, 'results/alignments/NC_018969_(.+).bam', '\\1'))
long_data

• In order to plot for a single isolate we need to first filter the data, then pass the data into ggplot(). Here we plot the plasmid position (POS) against the number of reads covering that position (depth) using geom_line() setting the line width to 2

long_data %>%
  filter(isolate == 'SEQ045') %>%
  ggplot(aes(x=POS, y=depth)) + 
    geom_line(size=2) 

• lets add two arrows to reprsent the position of the genes on the plasmid, we just prove x,y coords and some style attributes

long_data %>%
  filter(isolate == 'SEQ045') %>%
  ggplot(aes(x=POS, y=depth, group=isolate)) + 
    geom_line(size=2) +
    geom_segment(aes(x=1,y=1500,xend=477,yend=1500),size=3,color='blue',arrow=arrow(length=unit(0.4,"cm"))) +
    geom_segment(aes(x=1292,y=1500,xend=2026,yend=1500),size=3,color='red',arrow=arrow(length=unit(0.4,"cm"))) 

• Finally we can change the axes labels with labs() and set the theme to theme_bw()

long_data %>%
  filter(isolate == 'SEQ045') %>%
  ggplot(aes(x=POS, y=depth, group=isolate)) + 
    geom_line(size=2) +
    labs(x='LlaG1 gene position', y='Sequence depth') +
    geom_segment(aes(x=1,y=1500,xend=477,yend=1500),size=3,color='blue',arrow=arrow(length=unit(0.4,"cm"))) +
    geom_segment(aes(x=1292,y=1500,xend=2026,yend=1500),size=3,color='red',arrow=arrow(length=unit(0.4,"cm"))) +
    theme_bw()

• Note that it is sometimes useful to plot the y axis with a log scale, using the function scale_y_log10()

long_data %>%
  filter(isolate == 'SEQ045') %>%
  ggplot(aes(x=POS, y=depth, group=isolate)) + 
    geom_line(size=2) +
    scale_y_log10() +
    labs(x='LlaG1 gene position', y='Sequence depth') +
    geom_segment(aes(x=1,y=1500,xend=477,yend=1500),size=3,color='blue',arrow=arrow(length=unit(0.4,"cm"))) +
    geom_segment(aes(x=1292,y=1500,xend=2026,yend=1500),size=3,color='red',arrow=arrow(length=unit(0.4,"cm"))) +
    theme_bw()

---

Plotting multiple charts

• if we want to plot the same graph for all isolates we can use facet_wrap(). not the removal of the filter() function and the inclusion of the group=isolate in the ggplot() function, then we apply facet_wrap(~isolate) to complete the process

long_data %>%
  ggplot(aes(x=POS, y=depth, group=isolate)) + 
    geom_hline(yintercept=0, color="red", size=2) + 
    geom_line(size=2) +
#    scale_y_log10() +
    labs(x='LlaG1 gene position', y='Sequence depth') +
    geom_segment(aes(x=1,y=5000,xend=477,yend=5000),size=3,color='blue',arrow=arrow(length=unit(0.4,"cm"))) +
    geom_segment(aes(x=1292,y=5000,xend=2026,yend=5000),size=3,color='red',arrow=arrow(length=unit(0.4,"cm"))) +
    geom_hline(data = coverages.NC_017763, aes(yintercept = coverage), linetype="dashed", color="darkgreen", size=1) +
    theme_bw() +
    facet_wrap(~isolate)

• We can save that final figure to a file with ggsave()

ggsave('NC_018969-mapping.png', height=20, width=20)

---

Exercise

• load the dataset stored in data/data.csv, plot a bar chart of county against popadults, but plot for all states

• set the y-axes to use a log scale
• change the y-axes labels to use commas in the numbers
• add a regresson line
• colour the points by the category column
• include all the subplots in one row. Hint: check ?facet_wrap

---

Session details

• Generate document version details

sessionInfo()

R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS:   /opt/R/4.1.1/lib64/R/lib/libRblas.so
LAPACK: /opt/R/4.1.1/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8    LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggsignif_0.6.3  ggdendro_0.1.23 gapminder_0.3.0 gganimate_1.0.7 plotly_4.10.0   patchwork_1.1.1 DT_0.23        
 [8] readxl_1.3.1    forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7     purrr_0.3.4     readr_2.0.1     tidyr_1.1.3    
[15] tibble_3.1.4    ggplot2_3.3.5   tidyverse_1.3.1 RSQLite_2.2.14 

loaded via a namespace (and not attached):
 [1] nlme_3.1-152      fs_1.5.0          lubridate_1.7.10  bit64_4.0.5       progress_1.2.2    httr_1.4.2       
 [7] tools_4.1.1       backports_1.2.1   bslib_0.3.1       utf8_1.2.2        R6_2.5.1          DBI_1.1.1        
[13] lazyeval_0.2.2    mgcv_1.8-36       colorspace_2.0-2  withr_2.4.2       tidyselect_1.1.1  prettyunits_1.1.1
[19] bit_4.0.4         compiler_4.1.1    cli_3.0.1         rvest_1.0.1       xml2_1.3.2        labeling_0.4.2   
[25] sass_0.4.1        scales_1.1.1      digest_0.6.27     rmarkdown_2.11    pkgconfig_2.0.3   htmltools_0.5.2  
[31] dbplyr_2.1.1      fastmap_1.1.0     htmlwidgets_1.5.4 rlang_0.4.11      rstudioapi_0.13   jquerylib_0.1.4  
[37] farver_2.1.0      generics_0.1.0    jsonlite_1.7.2    crosstalk_1.2.0   vroom_1.5.5       magrittr_2.0.1   
[43] Matrix_1.3-4      Rcpp_1.0.7        munsell_0.5.0     fansi_0.5.0       lifecycle_1.0.0   stringi_1.7.4    
[49] yaml_2.2.1        MASS_7.3-54       plyr_1.8.7        grid_4.1.1        blob_1.2.2        parallel_4.1.1   
[55] crayon_1.4.1      lattice_0.20-44   haven_2.4.3       splines_4.1.1     hms_1.1.0         knitr_1.34       
[61] pillar_1.6.2      reprex_2.0.1      glue_1.4.2        evaluate_0.14     data.table_1.14.0 modelr_0.1.8     
[67] vctrs_0.3.8       tzdb_0.1.2        tweenr_1.0.2      cellranger_1.1.0  gtable_0.3.0      assertthat_0.2.1 
[73] cachem_1.0.6      xfun_0.26         broom_0.7.9       viridisLite_0.4.0 memoise_2.0.1     ellipsis_0.3.2